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Abstract 

The proper specification of boundary conditions at artificial boundaries for the 
simulation of time-dependent fluid flows has long been a matter of controversy. 
In this work we apply the general theory of asymptotic boundary conditions for 
dissipative waves, developed by the author in [8], to the design of simple, accurate 
conditions at a downstream boundary for incompressible flows. For Reynolds num- 
bers far enough below the critical value for linear stability, a scaling is introduced 
which greatly simplifies the construction of the asymptotic conditions. Numerical 
experiments with the nonlinear dynamics of vortical disturbances to plane Poiseuille 
flow are presented which illustrate the accuracy of our approach. The consequences 
of directly applying the scalings to the equations are also considered. 


1 Introduction 

Problems posed on unbounded domains are ubiquitous in theoretical and applied fluid 
mechanics. Examples include exterior domains, as arise in aerodynamics, and cylindrical 
domains, as are found in the study of internal flows. Numerical simulations require a 
finite computational mesh which is typically obtained by introducing an artificial bound- 
ary. The proper specification of boundary conditions there has long been a matter of 
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controversy. Their choice is often made on the basis of ad hoc reasoning and their effect 
on the solution not fully determined. 

In [8] the author develops a general theory of asymptotic boundary conditions for 
linear, dissipative wave propagation problems. The basic approach is to numerically 
identify dominant wave groups and construct linear approximations to the dispersion 
relation which are valid in the neighborhood of the dominant waves. These lead to simple, 
local boundary operators. Error estimates for the truncated problem are derived which 
establish convergence as the size of the computational domain is increased. The goal of 
this work is to apply the general theory to the linearized, incompressible Navier- Stokes 
equations and to test it in fully nonlinear simulations. 

In Section 2 we outline the construction of the asymptotic expansions and boundary 
conditions for the linearized Navier-Stokes equations. Simplified approximations to the 
full theory, valid for moderate to large Reynolds numbers, are also constructed. These 
lead to boundary conditions involving fixed constants and a parametric dependence on 
Reynolds number. The simplifications are based on a simple scaling argument. In Sec- 
tion 3 we carry out the numerical computation of the asymptotic boundary operators 
for the special case of linearizations about plane Poiseuille flow. Both the simplified 
and full theory are considered, the latter to validate the approximations leading to the 
former. Numerical experiments with the nonlinear dynamics of vortical disturbances 
are presented. These verify the accuracy of our approach even for large amplitude per- 
turbations. Comparisons with the solution on an extended domain are given to study 
the decay of the error with increasing domain size. Finally, in Section 4, we develop 
an alternate approach based on the direct introduction of the scalings used to simplify 
the asymptotic analysis. This leads to a simple set of boundary conditions previouslv 
considered by Naughton [13]. 

Boundary conditions for similar equations based on the analysis of linearizations about 
homogeneous base flows have been derived by Halpern [9] and Halpern and Schatzman 
[10]. Related techniques have been employed by Abarbanel and coworkers [1] for the 
compressible Navier-Stokes equations, by Fasel [7] for simulations of transition and by 
Ache [2] for steady flows. More recently, Danabasoglu and coworkers [5] have investigated 
the effect of outflow conditions on the simulation of spatially evolving instabilities. 
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2 Construction of asymptotic boundary conditions 
for the incompressible Navier-Stokes equations 


2.1 Asymptotic analysis 


We consider 
flow: 


the incompressible Navier-Stokes equations linearized about plane parallel 

fc + tWfe-t/'Mg « £(& + &) 


dV d 2 rp 

^ 3x 2 <9y 2 ’ 

(*.y) € (0, oo) x (y 0 ,yi), 
dib n 

¥> = -r- = 0, y = y 0 ,yi- 

oy 

u’(0. y, t) = go(y,t), w(0,y,t) = yi(y,0- 


( 1 ) 

( 2 ) 

(3) 

(4) 

(5) 


Here, rp is the stream function perturbation, u; is the vorticity perturbation and L (y) is 
the parallel base flow. For generality we have included eddy viscosity terms consistent 
with a simple mixing length model. (See, e.g., Stanisic [14].) Rei is, then, a Reynolds 
number based on the mean flow and the kinematic viscosity while Re t is a turbulence 
Revnolds number. An unbounded domain in y could also be considered. 

Equations (1-5) are meant to model the ‘far field for a more complicated and gen- 
eral class of problems. Our goal is to derive reasonably simple but accurate boundarv 
conditions at a downstream boundary, x = t. The general approach to the construction 
of asymptotic boundary conditions given in [8] may be directly applied in this case. To 
begin we consider the eigenvalue problem: 


sjj + A U{y)Cj - A U"{y)rj} 


rT, + 9$) 

+^(D(y)^)-£D"(y)rp, 


( 6 ) 


H > — X 2 tp + 


d 2 Tp 

ly 2 ' 


{-) 


The equations above are well known in the theory of hydrodynamic stability and are usu- 
ally referred to as the spatial Orr-Sommerfeld equations. Techniques for their numerical 
and asvmptotic analvsis have been extensively developed. (See, e.g. Drazin and Reid 

[ 6 ])- 
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( 3 ) 


Solutions to the linearized problem may be expressed in terms of pairs, 



satisfying 9fc(A/(s)) < 0 for 9?(s) sufficiently large. We have: 

( 5 )-?*■• 

= f ci(p)q,(x,y,t - p)dp, 

J 0 

= -T-. j t )t+XMx ( \ d s _ 

2m Jc \ ip t (y\s) J 

Here C is an appropriate inversion contour. 

By a standard steepest descent computation we formally obtain: 


( 9 ) 

( 10 ) 

(11) 


( w,(y;s-( 7 )) \ 

v /2tta;'(5-(7))^ V My;s m h)) J 


( 12 ) 


where 

t 

7 = ~, x > 1, (13) 

and s" satisfies: 

A;(s-) = -7, »(7)>0, 3(7) = 0 . (14) 

Substituting the asymptotic expansion of qi into (10) yields an approximation to $/. 
The expansion itself may be interpreted as a description of wave packets moving at their 
group velocity, A. The formula for is then a superposition of wave groups generated 
at various times. 

In contrast with the usual problems in dispersive wave propagation, our generalized 
wave groups decay (or grow) exponentially as they propagate with a rate: 


Ml) = £(*‘(7)7 + M s *(7))) # 0 . (15) 

We exploit this property to restrict our attention to a small number of dominant wave 
groups. That is, we expect the signal far downstream to be dominated by wave groups 
for which pi is maximized. Differentiating with respect to 7 and using (14) we find a 
necessary condition for a local maximum: 


Ml) = *(*’(7)) = 0 . 


(16) 
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Let 7 1 correspond to a global maximum of /*/. The exact, nonlocal boundary condition 
satisfied by qi is given in transform space by: 

( - A,w) (17) 

We approximate it by replacing A; by the linear part of its Taylor series at s (7/)- This 
leads to a local condition: 

(£ " M*‘(7l)) -Ai(<‘(7i))( jj - *1 = O' (18) 

This operator may be directly applied to ^ ^ ^ if it corresponds to a maximum over 

all /. More generally, an operator may be constructed as the product of operators taken 
from a few modes and local maxima. 

We summarize below a procedure for constructing the asymptotic boundary condi- 
tions: 

• Compute local minima of the decay rates, —pi- That is compute solutions of (14) 
and (16). 

• Choose from these a small subset, (5,-, l } ), including the overall minimum decay rate 
as well as any other wave groups whose decay rate is close enough to the minimum. 


Use as the boundary condition at x = r: 

A (fo ~ A b( 5 i) - ~ 



= B 





(19) 


In [8] an error analysis is developed for this procedure. The main result is that the 
difference, e, between the solution on the truncated domain and the true unbounded 
domain solution satisfies: 

||e|| = 0 (^~ ), t — + 00. ('-0) 

Here, 

Mo = max/x/(7). (21) 

Although this result was obtained under fairly general hypotheses, we have been un- 
able to verify them in this case. The error estimate is, however, consistent with the 

computational experiments we will present. 

In order to directly complete the search for the dominant wave groups, the repeated 
numerical solution of (6-7) in conjunction with a line search procedure is required. We 
have implemented an algorithm to accomplish this using a spectral discretization of the 
0 rr- Sommerfeld problem [4]. This will be discussed in greater detail in the context of 
the numerical experiments. 
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2.2 Simplified approximations 

Despite the availability of a number of appropriate techniques, the direct numerical ap- 
proach to the construction of the asymptotic boundary conditions has certain drawbacks. 
Among these are the expense of the computation and the need to recompute the boundary 
condition as the Reynolds number(s) change. Consequently, we have developed an alter- 
nate procedure based on the asymptotic analysis of (6-7) for Re t y large. The assumption 
underlying this analysis is: 

Assumption 1 The dominant wave groups satisfy: 


1 = 0 . 


( 22 ) 


VVe note that this condition does not necessarily hold for all Reynolds numbers. In 
particular, near the critical Reynolds number from linear stability theory the dominant 
wave group corresponds to the critical Tollmien-Schlichting waves which, in general, do 
not have s' = 0. For the special case of Poiseuille flow, which we have used for numerical 
experimentation, we find that Assumption 1 holds for a wide range of subcritical values 
of Re/. 

By itself Assumption 1 allows the construction of the asymptotic boundary conditions 
using a single solution of the Orr-Sommerfeld equations. Furthermore, the necessary 
condition for its validity, (14), can be directly verified. Further simplification occurs 
when we consider an expansion of the solution of (6-7) in powers of (For laminar 
flows this is replaced by an expansion in ^). We take: 


A = — 1 

Re t ’ Re t 


To leading order the eigenvalue problem becomes: 

d 2 


1 ("wjp - * = 

where we have made the substitution 


and 


u> = 


Cf = 


,dy 2 

d 2 ifi 

W' 

Re t 

Re\ 


d 2 d 2 ' 


dy 2 


0, 


(23) 


(24) 


(25) 


(26) 


We may now proceed to construct the asymptotic boundary conditons based on a 
numerical analysis of (24) with Assumption 1. We emphasize the advantages of this 
simplification: 
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• We need only solve (24) once for 5 = 0 and then compute the required derivatives, 

d\ 
d f 

• The resulting boundary conditions depend algebraically on the Reynolds number. 

Unfortunately, (24) with s = 0 is not self-adjoint, so it is difficult to obtain general 
results on its solutions. An ideal result would guarantee the existence of solutions with 
real, positive group velocity. We note that a real group velocity is guaranteed for real 
eigenvalues. It is tempting to rewrite the troublesome term, 35 'ts symmetric part 

plus a perturbation, but the unbounded ‘perturbation’ is too large for the direct appli- 
cation of the general theory (Kato [12]). It is possible to derive some simple equations 
for the real and imaginary parts of A: 


where 


»(A) = t jrr > 

{h + + £) 

(27) 

0(A) = — S(A )£, 
1 1 

(28) 

= [ (D(y) + ei)\xl>"\ 2 dy, 

(29) 

h = r 

Jyo 

(30) 

h = \ r U"(yM\ 7 dy , 
*'Vo 

(31) 

3 = 9( r U'{y)Vi'dy). 

(32) 


These relations suggest that all but finitely many eigenvalues will have negative real part 
and that they cannot lie too far from the real axis. For the special case of Poiseuille 
fl ow . y 0 = -1, t/j = 1 and U = 1 — y 2 , our numerical solution of (24) indicates that 
all eigenvalues are real and negative and that the group velocities are positive. Again, 
these conditions may be checked in the general case by a single solution of (24). The 
computational cost of this procedure is then negligible in comparison with that of the 
full simulations. 


3 Numerical experiments: the dynamics of some 
vortical disturbances to plane Poiseuille flow 

In this section we describe the results of numerical experiments with the asymptotic 
boundary conditions we have developed. We take our base flow to be laminar channel 
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flow: 


U(y) = i - y 2 , y € (- 1 , 1 ), (33) 

and consider (1) without the eddy viscosity terms. (That is, D = 0 and Re = Rei.) 
We divide our discussion into two parts considering first the eigenvalue problems used to 
construct the asymptotic boundary conditions and, then, the use of these conditions in 
the simulation of the full nonlinear equations. 

3.1 Numerical analysis of the full and reduced eigenvalue prob- 
lems 

We have developed two separate programs for the computation of the various quantities 
defining the boundary conditions. The first solves a discrete approximation to the reduced 
eigenvalue problem (24) with s = 0 to find A and It is the simplified condition defined 
by this reduced problem which is used in the simulations. The second solves a discrete 
(pseudospectral) approximation to the full spatial Orr-Sommerfeld equations, (6-7), and 
also computes the group velocities. The results of the latter are used to validate the use 
of the reduced equations for moderate Reynolds numbers. In practice, the point of the 
approximations leading to (24) is to avoid the expense of solving (6-7), and we show that 
this is justified in the case of Poiseuille flow for Re large enough compared to 1 but far 
enough below Re crit v 5775. For larger values of Re the full eigenvalue problem must be 
used. 

Approximating (24) we introduce a uniform mesh, y, = — 1 + ih, i = 1 , . . . , A r , h = 
Denote by 't, the approximation to 0(y,). Using second order centered difference 
approximations to the derivatives we have: 


A (2tf, + (1 - yf)D + DA,) = (D + D_) 2 * t , 

(34) 

where 

(35) 

(D D 

/ i 4 

(36) 

These are supplemented by the discrete wall conditions: 


'f'o = 4v +1 = 0, 

(37) 

= #i, «iv +2 = 4 v. 

(38) 
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A 

dA 

tit 

-21.659 

-1.999 

-28.111 

-1.384 

-72.694 

-1.991 

-85.178 

-1.430 


Table 1: Dominant eigenvalues 


The eigenvalues and eigenvectors of the resulting generalized matrix eigenvalue prob- 
lem were computed using the IMSL subroutine EIGZF. To compute ^ we differ- 

entiate (24) with respect to s to obtain: 




dA, 

ds 


Lj>i- 


dy 2 ' 


(39) 


Here L and M are, respectively, the partial differential operators appearing on the left and 
right hand sides of (24) with s = 0. Multiplying this equation by an adjoint eigenfunction 
and integrating leads to an expression for For the discrete approximation we simply 
use the original difference approximation to L and Z)+Z)_ for Left eigenvectors, 

T, were computed by inverting the matrix of right eigenvectors and were then used to 
compute 

— Lili£±£zli — (40) 

ds f , (2$, + (1 - yfiD+D-Vi) 

All eigenvalues were found to be real and negative. Tabulated are the six smallest eigen- 
values along with the corresponding values of For these we have N = 39. Their 
accuracy may be verified by comparing with solutions of the full Orr-Sommerfeld equa- 
tion for s — 0 published by Bramley and Dennis [3]. 

The first and third eigenvalues correspond to even stream function perturbations 
while the second and fourth correspond to odd perturbations. It is interesting to note 
that the group velocities, -(^) _1 , are nearly equal within each symmetry class. In the 
experiments which follow a product boundary condition based on the first two modes is 
employed. This is formally justified if e V. > e~ However, the near equality of 
the group velocities is likely to extend the range for which the two mode condition is 
reasonably accurate. 

It is of some interest to view the flow field which results from adding to Poiseuille 
flow the asymptotic solution corresponding to the first symmetric and antisymmetric 
eigenfunctions. These are shown in Figure 1 for Re = 100, c = A6(t — to), t = to + 3. 
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It must be emphasized that the amplitude of the perturbation, A, is here chosen to be 
much larger than a linearized analysis can reasonably allow. Nonetheless, the streamlines 
bear interesting similarities to those we have obtained in full nonlinear simulations. 

Although the solutions of the reduced eigenvalue problem described above satisfy 
necessary, local conditions for dominant wave groups, sufficient conditions are global 
and more difficult to check. Indeed, it is well known that Poiseuille flow becomes linearly 
unstable for Re « 5775 [6]. Then the dominant wave groups are expected to grow in 
space, while those above decay. In order to further investigate the dominance of the 
solutions of the reduced problem for subcritical Reynolds numbers, as well as compute 
the asymptotic boundary conditions for higher Reynolds numbers, approximate solutions 
of the full spatial Orr-Sommerfeld equation were computed. 

The mesh points are taken to be the Chebyshev nodes: 

y, =cos( ^ - ~ 1) J, i = l,...,N + l. (41) 

The ijth component of the (.V — 1) x (;V — 1) differentiation matrix D\ is defined by: 

• Let Pj{y) be the polynomial of degree N which is 1 at y = y J+ \ and 0 at the other 
nodes. 


• Let (D i), 7 = Pj(y, + i), ij = 1 , . . . , JV - 1 . 

The zjth component of the (N + 1) x (JV + 1) differentiation matrix D 2 is defined by: 

• Let Pj{y) be the polynomial of degree .V which is 1 at y = y : and 0 at the other 
nodes. 


. Let (D 2 ), ] = P'/(y l ),i.j = 1,....V+1. 

Take U2 to be an (;V -f 1)- vector whose zth component approximates c j(y,) and = 

A u 2 . Let U3 be an (N — l)-vector whose ith component approximates At^(y 1+1 ) and u 4 
the (N — l)-vector whose zth component approximates ip'(y,+ 1). Note we are explicitly 
imposing the boundary conditions that ip and its normal derivatives be zero at y = ±1. 
A discrete approximation to the eigenvalue problem (6-7) is then given by: 


/ Re ■ V — D2 + Re • s • / Re ■ T 0 \ 


V 


I 

0 

0 


0 

w 

0 


0 

0 


Here the (N + 1) x (N - 1) matrix T is given by: 


T = 



0 

D , 

0 / 


u = Xu. 


(42) 


(43) 
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Figure 1: 
The signa 












while the ( N — 1) x ( N + 1) matrix W is given by: 


w=(o i o). 


( 44 ) 


The problem was solved for A given s using the IMSL subroutine EIGCC. Approximations 
to dj were computed following the same ideas used for the simplified eigenproblem. 

A (coarse) sampling of Reynolds numbers between 100 and 6000 was taken. For each 
value of Re solutions for s = k = 1,...,21 were found. This choice is motivated 

by (16) and the fact that neutral modes in the classical theory satisfy |£$(s)| < 1 [6]. If 
the imaginary part of the derivative of any eigenvalue appeared to pass through 0 and 
the real part of A was of the same order as the solutions of the reduced problem, the s 
increments were refined in the appropriate neighborhood. Note that a real corresponds 
to a maximum or minimum of 3f(A) as s varies on the imaginary axis. On the basis of 
these computations we tentatively conclude that the modes approximated by the reduced 
eigenvalue problem are dominant for a range of Reynolds numbers including at least the 
interval 100 < Re < 4000. 

It is somewhat surprising, and perhaps of some interest in the study of subcritical 
transition phenomena, that the reduced problem still provides the dominant modes at 
so high a Reynolds number. For example, for Re = 4000 we find that the least damped 
mode with s' 0 has a decay rate more than one and a half times that of the least 
damped mode with s =0. For Re — 6000, on the other hand, there is an eigenvalue 
with positive real part for 3(s) in an interval including (.25,-28). The dominant wave 
group is then described bv 

s* « .27z, (45) 

A % 9.11 x 10~“ - 1.027i, (46) 



(47) 


3.2 Numerical method 

In order to test the proposed boundary conditions, we consider the following problem: 

d<jj drl’ du di> du \ ( d 2 u> d 2 u> \ 

lh ~ l^lh + Ihfy = Ye \d^ + d^)' 


_ d 2 ip <9 2 0 
{x,y) 6 (0, t) x (-1, 1), 


(49) 

(50) 
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with boundary and initial conditions: 


(51) 

(52) 

(53) 

(54) 

(55) 

(56) 

(57) 

y; — — 1 + ih, z = 1, . . . , JV; I, = jh, j = 1, . . . M, (58) 

was introduced and the spatial derivatives replaced by the standard second order central 
differencing formulas; D+D _ for the second derivatives and D 0 for the first. Time differ- 
encing was Crank-Nicholson for the viscous terms and Adams-Bashforth for the transport 
terms so that for interior nodes we have: 

^ _ _|_v 2 j J t+6t) =(l + - y(3 r (,) - T {t ~ St) ), (59) 

jt+st) = ( 60 ) 

Here V 2 is the discrete Laplacian and T ^ is given by: 

r (<) = -Doy^Dor^ + Dor^DoyJ'l (61) 

For all our simulations we have taken = .125. The wall boundary conditions on 0 are 
also directly discretized and u at the wall is computed by (60). The asymptotic boundary 
condition (54) was replaced by a product of discrete boundary operators following Higdon 
[11], A typical term, 

B,= 


a, 

dx Q dt ) ' 


0(x,y,O) = -y + y, u(z,y,0) = 2y, 


0(X,-1 ,t) = |. 0(X,1,O = -|» ^(x,±M) = °i 
0(0, y,t) = g 0 (y,t), w(0 ,y,<) = gi(y,t), 

50(r. y,t) = B ^-y + , Buj(t, y, t) = B(2y). 

Here the asymptotic boundary operator B is given by: 

»-(£->■-£ ■) (I -*■-*)• 

Ai = _ 2I | — , ai = -1.9990265, 

Re 

A 2 = ~ 28,1113 ^ , a 2 = -1.383905. 

Re 

A uniform mesh, 
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is approximated by: 



u( v + + v +~ v -~ *C) ~ + v+ + v+ + v_) 

-Wt( v t + v -- v +~ v -)^ 


( 63 ) 


where 



V) 

<**♦)■ 


(64) 


Then n^J becomes l\Bi. In the present case we have the product of two first order 
operators. The resulting discrete operator then involves three mesh lines in x and three 
time levels. 

To solve the linear system required by the implicit time stepping we employed what 
we believe to be a novel direct method based on discrete separation of variables. That 
is a method similar to fast solvers for elliptic equations. We begin by diagonalizing a 
rearrangement of the 2N x 2 N matrix: 


B - ( W 0 \ 

l / ~D +y D.J- 


(65) 


Although B is not symmetric, we find its eigenvalues to be real and positive for 
sufficiently large, certainly within the range of problems we have considered. Denote its 
eigenvalues by with corresponding right eigenvectors f t . Writing, 



( 66 ) 


we obtain an uncoupled collection of tridiagonal systems for the c l , s: 

K - = /„ 


(67) 


where the /,’ s are computed from the right hand sides of the implicit time stepping 
formulas. Note that the boundary conditions at 0 and r are directly applied to c,. 
A formal operation count for this method results in an estimate of 8N 2 M + o(N 2 M) 
operations per time step, which is the the same as for back substitution with a band 
solver (without pivoting). However, the preprocessing is likely to be cheaper for our 
method, at least when M is large enough compared with N (long channels). 


3.3 Results 

We first present results for Re = 400. We take N = 39 and consider r = 2,4,6. 
(M = 40,80, 120.) For purposes of comparison, we also compute a solution with r = 15, 
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Figure 2: Vorticity contours at t = 2.5, r = 15. 


which we will refer to as ‘exact 1 , and for r = 4 with Neumann boundary conditions, 
=z = 0. The inflow perturbations are taken to be: 



Xe- 6 <* a+ ( 1 - 2, > 2 >, 0 < t < 1 
0, t > 1 


( 68 ) 


9i = 2 y + 


24.4e“ 6(y2+(1-2f)2) , 0 < t < 1 
0, t > 1 


(69) 


Here A is an adjustable amplitude parameter. 

Figure 2 depicts vorticity contours at t = 2.5 for r = 15 and A = ^. Note that the 
flow certainly cannot be considered a small perturbation of Poiseuille flow, so that we 
cannot expect the linearized analysis leading to the boundary conditions to be valid. We 
are nonetheless interested in testing the conditions in this nonlinear setting. Figures 3-5 
depict the same flow at t = 6 as well as simulations with r = 4 using the asymptotic 
boundary conditions and Neumann conditions. We see that the exact solution and the 
one computed with the boundary conditions we propose are virtually indistinguishable, 
while the flow computed with the Neumann conditions is greatly distorted. Note that the 
displacement of the vortex causes a large error at the upper wall which extends well back 
into the computational domain. Figures 6-8 depict the same computations at t = 7 .5. 
Although the main disturbance has passed outside the computational domain, the error 
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Figure 3: Vorticity contours at t = 6, r = 15. 
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Figure 4: Vorticity contours at t = 6, r = 4, asymptotic conditions. 
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Figure 5: Vorticity contours at t = 6, r = 4, Neumann conditions. 
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Figure 7: Vorticitv contours at t = 7.5, r — 4, asymptotic conditions. 
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Figure S: Vorticity contours at t = 7.5, r = 4, Neumann conditions. 
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A 

Boundary Condition 

r 

Maximum Error 

~r 

7 

Asymptotic 

2 

1.8 

— n 

7 

Asymptotic 

4 

.65 

7 

Neumann 

4 

7.0 

—r~ 

7 

Asymptotic 

6 

.42 

-f- 

Ifi 

Asymptotic 

2 

.12 

T 

1 fi 

Asymptotic 

T 

.018 

~X" 

is 

Asymptotic 

6 

.013 ”1 


Table 2: Vorticity errors for 0 < t < 10. 

at the walls in the Neumann computation is still quite evident. 

It is also interesting to look at the maximum errors as a function of r. Here, again, we 
take the error to be the difference between the given solution on the truncated domain 
and the solution for r = 15. The tabulated results clearly show the decay of the error 
with increasing r. Indeed the theory predicts a decrease in error of approximately | as 
t increases from 4 to 6. and the results show this precisely. This conclusion can only be 
tentative, however, as at the end of the simulation the disturbance had not fully passed 
through the boundary in the case of r = 6. The error decreases much more rapidly as r 
increases from 2 to 4, but this is not surprising as our entire construction is dependent 
on large r asymptotics. 

Nonlinear effects are demonstrated by comparing results for A = ^ with those for 
.4 = i Although the amplitude is scaled by a factor of S, the errors scale by factors 
of 15 to 30. Nonetheless, we are quite surprised by the success of our method in the 
nonlinear regime. It cannot be explained by our linear analysis. It seems a nonlinear 
analysis of disturbance propagation is required. This is an interesting topic for future 
research. 

A limited number of simulations have also been carried out at Re = 2000 for the same 
inflow disturbances. Figures 9-10 show vorticity contours at t = 5.75 for simulations 
using, respectively, the asymptotic and Neumann conditions at r = 4. In this case we 
have N = 79 and M = 160. Due to the expense involved, we were unable to carry out a 
similar simulation on a large domain. Nonetheless, the differences in the vortex dynamics 
for the two downstream conditions is clear and, apparently, analagous with the results 
at lower Reynolds number. 
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Figure 9: Vorticitv contours at t = 5.75, r = 4, asymptotic conditions. 
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Figure 10: Vorticitv contours at t = 5.75, t = 4, Neumann conditions. 
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4 Pri mi tive Variable Formulation and Direct Ap- 
plication of the Scalings 


In this section we discuss the application of the ideas in this paper to the derivation 
of asymptotic boundary conditions for the primitive variable formulation of the Navier- 
Stokes equations. We are guided by the theory of boundary conditions for primitive 
variable formulations which is developed in the doctoral dissertation of Naughton [13]. 
This theory is based on the construction of energy estimates for the linearized equations. 
We emphasize that the conditions we will describe have not yet been fully analyzed or 
implemented. 

Linearizing about a parallel flow and assuming, for turbulent flows, a simple eddy 
viscosity we have: 

p- + p- = °. (™> 

ox oy 


f“ + £%)£ + U'W* 




dy ' Rei \9r 2 &V 2 / 

dp j 1__ ( d 2 v i d 2 v \ 

dy * Rc\ \5x 2 dy 2 ) 

+£, {D(y)£k + i( D M §;))■ 


( 71 ) 

(72) 


It is not entirely clear from the form of the equations what the correct number of 
boundary conditions is at an outflow boundary. However, following either the analysis 
of Halpern and Schatzmann [10] or Naughton [13] we seek two conditions. Often in nu- 
merical algorithms the somewhat awkward divergence equation is replaced by a Poisson 
equation for the pressure. Then, an additional boundary condition is needed. To guar- 
antee that the velocity field is divergence free, this condition must in general ensure that 
the divergence is zero at the boundary. (For outflow boundaries Naughton has shown 
that this may be relaxed to a Neumann condition on the divergence.) 

One possibility is to apply the boundary operator derived in the sections above to 
two of the perturbed variables. For example, 


Bu = 0, 


(73) 


Bv = 0. (74) 

An alternate approach, which could also be applied to the stream function - vorticity 
formulation, is to introduce directly into the equations the scalings used to approximate 
the spatial Orr-Sommerfeld equation. In terms of the physical variables these are: 
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( 76 ) 


(Here Re is the turbulent or kinematic Reynolds number as appropriate.) The continuity 
equation, (70), may be integrated from a wall boundary to yield: 


Differentiating (77) with respect to x yields: 

dv 1 

while substituting it into (72) along with the scalings leads to: 

!-<.>■ 


(77) 


(78) 


(79) 


Dropping the O(^j-) terms results in equations which could be used as boundary 
conditions. Unfortunately, the pressure condition is difficult to analyze using energy 
methods. To rectify this we add to it another 0( ^y) term and integrate with respect to 

y- 

P “^ + RT t D ^lh + ( 8 °) 

This corresponds to the so called normal constraint of [10]. The function $(f) corresponds 
to the arbitrary constant which may be added to p. For simplicity we take $(<) = 0. 
This leads to the following boundary conditions: 



(81) 


,1 1 . ..du 

P ~ ^Rei + Rei D ^dx = °' 

For solvers which use the Poisson equation for p this must be supplemented by either 

S + 8 =■«<*£» -»• 

Equations (81-82) are discussed in [13] where they are called Type 3 conditions. 
Energy estimates independent of Re are given there for the half channel laminar flow 
problem ( D = 0) using these conditions and linearizations about an arbitrary velocity 
field with outflow. A simple generalization of Naughton’s estimates holds in the case of 
turbulent flow. 

It is interesting to compare (81-82) with the conditions given in [10]. The equation 
for v is common to both works. In [10] equation (82) is replaced by a more complicated 
relationship between p and the velocities which is nonlocal in space. 
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In the future we plan to carry out numerical experiments which test these conditions 
as well as to consider the problem of deriving error estimates. Among other interest- 
ing questions to be studied are the application of these ideas to viscous, compressible 
flows, their use in the simulation of interesting unsteady flows and their generalization to 
problems in three dimensions. We note that the reduction of domain size allowed by the 
use of asymptotic boundary conditions is of increasing importance for three dimensional 
simulations. In the three dimensional case the reduction in effort resulting from the use 
of the approximate eigenvalue problem or the direct use of the scalings would also be 

most keenly felt. 
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